#!/usr/bin/perl
#########################################################
# folding.pl                                            #
#                                                       #
# use:    folds a time series to find periodicities     #
#                                                       #
# input:  a time series contained in a text file which  #
#         is arranged in colums such that the fist      #
#         column is a list of times and the second      #
#         column is a list a list of amplitudes         #
#         correspoding to thoes                         #
#         times                                         #
#                                                       #
# output:                                               #
#                                                       #
#########################################################
# code adapted from that writen by KMB                  #
# BM 10/12/04                                           #
#########################################################

#use warnings;
#use strict;

use Math::BigFloat;
use POSIX;
use Math::Trig;

BEGIN {
#========== varables ===================================

    use constant PI => 3.141592653589792384;
    $rcount = 0;         # number of data points
#   $frequency = 31.0;  # fold frequency
#   $frequency = 10.0;
#   $frequency = 0.17;
#   $p = 1/$frequency; 
#   $p = 13.08;
#   $p = 13.08211;
    $p = 0.04;
    
#   print "frequency: $frequency\n";
    printf("period: %5.3f\n",$p);
    $epoch  = 2450023.62;         # time origin offset
    $epoch = 2450023.62 - 485.*$p;
    $epoch = 0;
    $phase  = 0;         # phase of the signal
                         # relative to the origin of the
                         # of time series / (2*pi)
    $j      = 0;         # indexing variable

    
    $bin_step = 0.1;    # bin width in phase units
                        # 0.05 gives 20 bins

#========= read data ==================================
    
  $infile = "fold_test.dat";
    open(INFILE,"$infile") || die "cannot open: $infile\n\n";

    while(<INFILE>) {
	split;
 	$time[$j] = $_[0];
	$amp[$j] = $_[1];
  	$phase[$j] = (($time[$j]-$epoch)/$p - int(($time[$j]-$epoch)/$p));
 	$readphase[$j] = $_[6];
	$j++;
    }
    close(INFILE);
    $rcount = $j; 

     print "\ndata read in\n";
     print "----------------------------------\n";
     print "number of lines of data $rcount\n";
     print "----------------------------------\n";
   
    for ($j=0;$j<$rcount; $j++) { 
   	    printf(STDERR "time: %.3f\tamp: %.3f\tphase: %.3f\treadphase: %.3f\n",
   	       $time[$j],$amp[$j],$phase[$j],$readphase[$j]);
    }    
    
#========== sum data into bins ===========================
    $NumberOfBins = 1/$bin_step;
    $ibin = 1;     # bin index starts at 1

    for ($step = 0; $step < 1; $step += $bin_step) {
    	for ($j = 0; $j < $rcount; $j++) { 
#    	    printf(STDERR "Bin %d %f -> %f\n", $ibin, $step, $step + $bin_step);
	        if ($phase[$j] > ($step + 0.0) && $phase[$j] <= ($step+$bin_step)) {
    	    	$bintot[$ibin] += $amp[$j];
	        	$countbin[$ibin]++;
	    }
	}
#	print "bin ($ibin) count: $countbin[$ibin] total: $bintot[$ibin]\n";
	$ibin++   # increment bin index - node the final $ibin is 1
                  # greater than the last index
    }

#========= caluclate mean amplitude in each bin ==============
    for ($j=1; $j<=$NumberOfBins;$j++) { 
	if ($countbin[$j] > 0.0) {
	    $binmean[$j] = $bintot[$j]/$countbin[$j];
	}
	printf (STDOUT "\$binmean[%2d] = %7.4f\n",$j,$binmean[$j]);
    }

#========= calculate the sumsqdev of each bin =================
    for ($step=0;$step<1;$step+=$bin_step) {
    	for ($j=0;$j<$rcount;$j++) { 
	        if ($phase[$j] > ($step+0.) && $phase[$j] <= ($step+$bin_step)) {
	    	    $sumsqdev[$j] += ($amp[$j] - $binmean[$j])**2.;
	    	    printf(STDERR "amp %f mean %f\n",$amp[$j], $binmean[$j]);
	        }
	    }
	$ibin++ # increment bin index - node the final $ibin is 1
	        # greater than the last index
    }

#========= calculate the standard deviation of each bin ======
    for ($j=1;$j<= $NumberOfBins;$j++) { 
	if ($countbin[$j] > 0.) {
	    $stddev[$j] = sqrt($sumsqdev[$j])/$countbin[$j];
	}
    }

#========== print output to STDOUT ===========================   
#   print "bin |  mid phase | bin mean | bin stdev\n";
#   print "-----------------------------------------\n";
    for ($j=1; $j<=$NumberOfBins;$j++) { 
	printf (STDOUT "%2d     %7.4f    %7.4f     %7.5f  %d\n",
		$j,(($j*$bin_step)-($bin_step/2.0)), $binmean[$j],$stddev[$j], $countbin[$j]);
    }    
}

